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The Kitaev-Hubbard model of interacting fermions is defined on the honeycomb lattice and, at 
strong coupling, interpolates between the Heisenberg model and the Kitaev model. It is basically 
a Hubbard model with ordinary hopping t and spin-dependent hopping t'. We study this model 
in the weak to intermediate coupling regime, at half-filling, using the Cellular Dynamical Impurity 
Approximation (CDIA), an approach related to Dynamical Mean Field Theory but based on Pot- 
thoff’s variational principle. We identify four phases in the (U, i') plane: two semi-metallic phases 
with different numbers of Dirac points, an antiferromagnetic insulator, and an algebraic spin liquid. 
The last two are separated by a first-order transition. These four phases all meet at a single point 
and could be realized in cold atom systems. 


I. INTRODUCTION 


Mott insulators are systems that should be metals 
within band theory, but are in fact insulators because 
of electron-electron interactions. 1,2 However, the Mott 
phase is often hidden behind a magnetically ordered 
phase at low temperature. 3 Spin liquids are non-magnetic 
Mott-insulators, without broken lattice symmetry, stabi¬ 
lized purely by quantum effects. 4 In addition to a spectral 
gap, they are characterized by spin correlations that de¬ 
cay either exponentially, or as a power law, in the case 
of algebraic spin liquids. 5,6 Experimentally, a spin liq¬ 
uid ground state has been suggested in the organic ma¬ 
terial k-(BEDT-TTF) 2 Cu 2 (CN) 3 , 7 in other systems like 
YMnCb 8 and, more recently, in materials with a kagome 
lattice structure. 9,10 Theoretically, spin liquid phases 
were found, for instance, in the spin-f Heisenberg model 
on the kagome lattice, 11-13 and in the intermediate¬ 
coupling Hubbard model on a triangular lattice. 14 

Tikhonov et al. 15 have shown that an algebraic spin 
liquid is realized when a special type of perturbation is 
added to the Kitaev spin model. 16 It can be shown that 
the stability of the spin liquid phase, in that system, is 
due to time-reversal symmetry. The existence of an alge¬ 
braic spin liquid in a model of interacting fermions, the 
Kitaev-Hubbard model, was shown by Hassan et al. 1[ 
The phase diagram of this model was investigated us¬ 
ing the variational cluster approximation (VCA) which 
allowed the authors to identify a semi-metallic phase, a 
Neel phase and an algebraic spin liquid phase. 

In this work, we refine the analysis of Ref. 17 by using 
the cluster dynamical impurity approximation (CDIA). 
This method is more accurate in its treatment of the 
Mott transition, which appears clearly as a discontinu¬ 
ous transition with hysteresis. We also reveal a topologi¬ 
cal transition within the semi-metallic region, between a 
phase with eight distinct Dirac points, and another phase 
with only two Dirac points; this is a Lifshitz transition, 
that carries into the interacting region. These four phases 


(the algebraic spin liquid, the antiferromagnet, and the 
two semi-metallic phases) meet at a single point in the 
phase diagram. 

The paper is organized as follows. We review the model 
and describe its non-interacting solution in Section II. In 
Section III, we review the methods used in the interacting 
case (the CDIA), before presenting and discussing our 
results in Section IV. 


II. THE NON-INTERACTING LIMIT 


We will mostly follow the notation of Ref. 18. The 
Kitaev-Hubbard model, defined on the honeycomb lat¬ 
tice, has the following Hamiltonian: 

H = E j c ! ) c > +H.c.|+U^n it n 4 (1) 


where Ci„ annihilates a fermion of spin a at site i (the spin 
index is implicit in the above), <r“ are the Pauli matrices 
(a = x,y,z), U the Coulomb repulsion for two electrons 
of opposite spin on the same site, ni a = c\ a Ci a is the 
number of electrons of spin a at site i , and (i,j) a denotes 
the nearest-neighbor pairs in the three hopping directions 
of the cluster system (see Fig. 4 below). Throughout this 
work we will express energies relative to t , i.e., we will 
set t = 1. 

In the limit where U t,t', and at half-filling, the 

Hamiltonian (1), becomes equivalent to a combination of 
the Heisenberg and Kitaev 16 Hamiltonians: 


H(2) = E 

(ij)c 
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( 2 ) 


Time reversal is applied by changing the sign of the Pauli 
matrices (<r“ -A — a a ). The Hamiltonian (1) breaks time- 
reversal symmetry explicitly when f / 0, but is invariant 
under parity (as defined by the exchange of sublattices 
A and B on the honeycomb lattice). 
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In the non-interacting limit (U = 0), the Hamilto¬ 
nian (1), after Fourier transform, can be expressed in 
terms of the destruction operators c^a and c^ b on the 
A and B sublattices (the spin index is, again, implicit): 


#0 — ^ ( c k,yt C k,s) 

k 





( 3 ) 


as t' > 0, a total of six new distinct Dirac points appear, 
along the lines that join K and K', i.e., on the Brillouin 
zone boundary (Fig. 1). These Dirac points form the 
Fermi surface at half-filling. 

The positions of the Dirac points can be found by solv¬ 
ing the equation £-(k x = k y ) = 0 for k y . In terms of 

A = cos(^ky), this is a quadratic equation: 


where S(k) = P 3 +P 1 e llC2 +P 2 e~ lkl is a 2x2 matrix acting 
in spin space, with the projectors P a = \( 1 + t'a a ) and 
&i( 2 ) = k • e 1( - 2 ). The vectors are a Bravais basis of 

the honeycomb lattice: e 1 ( 2 ) = (±§, ^). 

For a given k, the four eigenvalues of Hq are ±e±(k) 
with £±(k) = \ [£(k) ± t'|B(k)|]. We have defined 

3 

£(k) = -(1 + if 2 ) + cosfci + cos k 2 + cos k 3 (4) 

and the vector B with components 

Bi (k) = 1 — tf sin k\ + cos k 2 + cos k 3 
B 2 (k) = 1 + cos k\ — t' sin k 2 + cos k 3 
B 3 (k) = 1 + cos ki + cos k 2 — if sin k 3 (5) 

with k 3 = —k\ — k 2 . The four-component eigenvector 
associated with eigenvalue p£ p >{ k) (jf = ± and p = ±) 
will be denoted <F pp . 


A. Dirac points 
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FIG. 1. Position and chirality of the Dirac points at U = 0: 
open (blue) and filled (red) circles are the center of positive 
and negative circulations of Vy, respectively. As t' increases 
from zero, the new Dirac points drift in the direction indicated 
by the arrows, until they annihilate at the points marked 
by crosses, at the critical value t' c , while the graphene Dirac 
points stay fixed. The hexagon is the Brillouin zone. 


At if = 0 (the graphene limit), the Fermi surface con¬ 
sists of two distinct Dirac points K = (27 t/3, 27t/ 3\/3) 
and K' = (27r/3, — 27 t/ 3\/3) that are the focal points of 
Dirac cones at positive and negative energies. As soon 


4[A 2 (1 — t' 2 ) + A] + 1 + t' 2 = 0 (6) 


whose solutions are 

Ai = - ^ and A 2 


1 + t' 2 

2(1 — t' 2 ) 


( 7 ) 


Ai corresponds to the graphene Dirac point at the zone 
corner, whereas A 2 is an additional Dirac point for a given 
if; the other six Dirac points can be deduced from Eq. (7) 
by lattice symmetries. The limiting case A 2 = — 1, for 
t! < 1, corresponds to a critical spin-dependent hopping 
t' c = \/y/Z. At that value of if the six additional Dirac 
points merge pairwise and disappear at the midpoints be¬ 
tween zone corners, as illustrated on Fig. 1. This merging 
is discussed in more detail below. 

Note that the Dirac points are protected by parity. 
Adding a parity-non-conserving term, such as a staggered 
magnetization M, would change the Hamiltonian (3) into 


v (V C f ) ( ¥ 3(k ^ 

\ c k,A c k ,b) ^st(k) —M) 



( 8 ) 


and the corresponding energies would become £±(k) = 
\ [£(k) + M 2 ± t'|B(k)|], which never vanishes. Thus all 
Dirac points disappear if M ^ 0. 


B. The Pancharatnam-Berry curvature 


The Dirac points are singular points of the 
Pancharatnam-Berry (PB) curvature. The latter is given 

by 


ijw'(k) = ^0„$w'(k)t3„$w'(k) (q) 

47H 

where /r, v = 1,2 are two orthogonal directions in the 
Brillouin zone, e pl/ is the two-dimensional Levi-Civita 
tensor, and the eigenstates of the Hamiltonian (3) are 


with 


$pp'(k) 
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C'(k) 


( B 3 (k)T|B(k)| \ 
Bi(k) + iB 2 (k) J 


( 10 ) 


( 11 ) 


and C(k) is a normalization factor. The expression for 
the phase x(k) is known analytically for all values of t' 
but is to complex to reproduce here. 



3 


It can be shown that 

i^'(k) =p' 1 - [6(k) +6(-k)] + l -e^d^d vX ^) (12) 
where 

6(k) = —B(k) • a M B(k) x 0 v B(k) (13) 

o7T 

The first term of (12) is everywhere regular, but the last 
term is singular where the phase \ is ill-defined, which 
occurs when two bands cross at a given wave vector, i.e., 
at Dirac points. 



k x -K x 

FIG. 2. Circulation of Vy(k) around the Dirac point K and 
the new, satellite Dirac points, at t' = 0.4. 



FIG. 3. (color online). Profile of the dispersion around |G = 
(K + K')/2, one of the merging locations of the new Dirac 
points at t' = t' c . The dispersion is shown along the x (left) 
and y (right) directions. The various curves correspond to an 
array of values of t' close to, and around t' c . Precisely at t' c , 
the dispersion is linear in k x and quadratic in k y . For t' ^ t' c ’ 
the two Dirac points stand away from |G. 


The integral of the PB curvature over occupied wave 
vectors is the Chern number. At half-filling, the two 
lowest-energy bands contribute 1 and —1, respectively, 
to the Chern number, coming from the first term of 
(12). Each Dirac point will in addition contribute 


to the Chern number, coming from the second term of 
(12), i.e., from the circulation of V\/47r around it. Fig¬ 
ure 2 shows this circulation around the old and new Dirac 
points in the vicinity of K. From that figure, it appears 
clearly that the graphene Dirac point K will contribute 
whereas the new Dirac points have the opposite con¬ 
tribution; but this picture is reversed when looking at the 
Dirac points surrounding K 7 . Thus the Dirac point con¬ 
tributions sum up to zero over the whole Brillouin zone, 
since they occur in pairs with opposite chiralities. 

At half-filling, there is particle-hole symmetry and the 
explicit time-reversal breaking in the Hamiltonian does 
not induce any topology, because of zero total Chern 
number. We can think of this as an accidental restoration 
of time-reversal symmetry. 

Note that the t' = 0 line (the graphene limit) is sin¬ 
gular in this respect. At t' = 0, the bottom two bands 
become degenerate, and likewise for the top two bands. 
Thus the graphene Dirac points K and K 7 arise for each 
of the two spin bands, and the two spins make exactly 
opposite contributions to the PB curvature. The latter 
is thus identically zero everywhere, whereas for t' > 0 
the PB curvature is not zero, but its integral over the 
Brillouin zone (the Chern number) is. 


C. Lifshitz transition 

At the critical value t' cl the new Dirac points with op¬ 
posite chiralities annihilate pairwise, at wave vectors 
lying midway between zone corners. These merging wave 
vectors, indicated by green crosses on Fig. 1, are time- 
reversal invariant, since — is equivalent to ^G be¬ 
cause G is a reciprocal lattice vector. This topological 
phase transition cannot be described by the total Chern 
number, which does not change here. It is a Lifshitz 
transition, akin to what has been described in Refs 19— 
21. Precisely at this transition, the dispersion around 
the merged Dirac points is linear in one direction and 
quadratic in the other, a behavior qualified as semi-Dirac 
in Ref. 20 and illustrated on Fig. 3. Experimentally, this 
transition could be probed by changes in the tunneling 
probability between the valence and conduction bands 
during Bloch-Zener oscillations. 19 


III. INTERACTIONS 

A. The cluster dynamical impurity approximation 

In the interacting case, the Hamiltonian (1) must be 
treated within some approximation method. In this 
work, we use the cluster dynamical impurity approx¬ 
imation (CDIA), 22,23 closely related to the variational 
cluster approximation (VCA) and to the cellular dynam¬ 
ical mean field theory (CDMFT), 24 but more accurate 
in its rendering of the Mott transition. These methods 
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can be understood in the framework of Potthoff’s self¬ 
energy functional approach (SFA). 25,26 In this approach, 
the physical self-energy of the system is obtained via a 
dynamical variational principle, expressed by the Euler 
equation 


= 0 (14) 

The self-energy functional f2[S] is defined as follows: 

fi[£] = F[£] + Tr ln[-(Go 1 - £) _1 ] (15) 

where Go is the Green function of the non-interacting 
part of Hamiltonian (1), and F[X] = $[G[X]] — 
Tr(£G[£]) is the Legendre transform of the Luttinger- 
Ward functional <F[G], 2 ' defined as a functional of the 
Green function G. We use a matrix notation for the 
Green function and self-energy to emphasize that they 
are to be considered as matrices in the space of fre¬ 
quencies and degrees of freedom (e.g. sites and spin). 
The symbol Tr means a functional trace, i.e a sum over 
bands, wave vectors and frequency. At the stationary 
point of H[S], the value of the functional H coincides 
with the thermodynamic grand potential E — fiN. 



grand potential fl'[S(h)] of these solutions, and thus the 
full Potthoff functional can be expressed as: 

n[E(h)} = n'[E(/i)] 

+ Tr InHGo 1 - S(h))" 1 ] - Tr ln(-G \h)) (16) 

where G' is the known physical Green function of the 
reference system. This relation provides us with an exact 
value of the functional f2[X(/i)], albeit on a restricted 
space of self-energies X(/i') which are the physical self¬ 
energies of the reference Hamiltonian H'. If we introduce 
the notation V(w) = G^ 1 — Gq - 1 , we can rewrite the 
above as 

fi[£(/i)] = ft'[£(/i)] + Tr ln[l - V(w)G'(w)] (17) 

Generally, the reference Hamiltonian H' is based on 
a finite, periodically repeated cluster or set of clusters. 
This periodicity defines a superlattice, and a correspond¬ 
ing reduced Brillouin zone, smaller than the original lat¬ 
tice’s Brillouin zone. The reference Green function G' is 
then momentum independent, and the matrix V depends 
on a single momentum (i.e., is diagonal in momentum in¬ 
dices). Eq. (17) then reduces to the more explicit form 



FIG. 4. (color online). The two-cluster system used in this 
work as a unit cell. Shaded numbered circles are lattice sites 
and small numbered circles represent the bath orbitals. The 
second cluster is a spatial inversion of the first. 


Unfortunately, the Luttinger-Ward functional H> [G], 
and consequently its Legendre transform E[S], is not 
known explicitly. This leads to some approximations, in 
the weak-coupling regime, where $[G] is represented by 
a truncated sum of diagrams. The Hartree-Fock approx¬ 
imation is an example of such a truncation. The basic 
idea behind the SFA is that the functional F[E] is an 
universal functional of the self-energy. 25,28 This means 
that the functional form of F[£] is the same for a ref¬ 
erence Hamiltonian H' with the same interaction as H , 
but a different non interacting part. Typically, H' will be 
a small system, e.g., a cluster, whose solution is known 
numerically. Given the physical self-energy S(/i) for a 
family of reference Hamiltonians H' parametrized by h, 
the value of F[X(/i)] can be extracted from the known 


Sl{h) = n\h) + lndet |l - V (k< w)G'(w)] (18) 

k 

where now the matrices are ‘small’, i.e., their order is the 
number of degrees of freedom in the repeated unit and 
N is the (potentially large) number of lattice sites. 

In the VCA, variational fields are added within the 
clusters, in order to give room to possible broken symme¬ 
tries. By contrast, CDMFT does not add extra terms to 
the cluster, but instead uses a set of non-interacting, ficti¬ 
tious orbitals (the bath) that are hybridized to the clus¬ 
ter and represent its immediate physical environment. 
In CDMFT, the Potthoff functional (16) is not calcu¬ 
lated, and the solution is found instead by imposing a 
self-consistency relation between the cluster Green func¬ 
tion G' and the projection of the lattice Green function 
G onto the cluster. In CDIA, a bath system is intro¬ 
duced, just like in CDMFT, but the solution is found by 
solving the Euler equation (14), like in VCA. This allows 
us to also introduce variational fields on the cluster if 
needed, and at the same time gives a better description 
of temporal fluctuations (because of the presence of a 
bath), which is important to correctly capture the Mott 
transition. 


B. Reference system 

The reference system used in this work is based on a 
unit cell made of two four-site clusters, the second ob¬ 
tained from the first by a spatial inversion and a shift, as 
illustrated on Fig. 4. Together, these two clusters form 
an 8-site supercluster that tiles the honeycomb lattice. 
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FIG. 5. (color online). Phase diagram of the half-filled 
Kitaev-Hubbard model (Eq. (1)) on the U — t' plane (t = 1). 
n_o is to the number of Dirac points in the semi-metallic (SM) 
phase; the t' = 0 limit corresponds to graphene (no = 2), rep¬ 
resented by a green line. The dashed line indicates the first- 
order Mott transition in the normal (non magnetic) state. 
The antiferromagnetic (AF) insulator phase is bounded by 
the red curve. The algebraic spin liquid (ALS) phase is the 
region of gapped spectrum and zero staggered magnetization. 



U/t 


FIG. 6. (color online) The antiferromagnetic order parameter 
as a function of interaction U for different spin-dependent 
hopping t'. We observe a continuous transition for 0 < t' < 
0.4, a weakly discontinuous transition when 0.4 < t! < 0.83 
and a discontinuous transition for t' > 0.83. 


Each of these two clusters contains four spatial sites and 
6 bath sites. The bath sites have no position per se, 
but can be thought of as representing the nearest sites of 
each cluster’s environment, and are illustrated with this 
in mind on Fig. 4 (small colored circles). This is why 
they are hybridized with the cluster boundary sites only, 


not the central site. 

The reference Hamiltonian has the following expres¬ 
sion: 


H ' = E | c i ( ^ + 2 g ) c i + Hc } + ^E n A n A 

+ E e /J, a \i<7 a n<j + E C i + H - C - 

fi,a 

(19) 


where a^ a annihilates an electron of spin a at the bath 
orbital /i, is the energy of the bath orbital /z, 0^,, is the 
hybridization between the bath orbital n and site z, and 
•di^ a corresponding spin-dependent hybridization. The 
Pauli matrix appearing in the bath hybridization 

is determined by the corresponding orientation (x , y or 
z) of the hybridization link on Fig. 4. e M , 9i M and i9^ will 
be treated like variational parameters, and the solution 
adopted will be such that Q(e,9) is stationary. Here the 
sum over i,j is restricted to the cluster. At a particle- 
hole symmetric point, such as the normal phase at /z = 
f//2, only two of these bath parameters are independent, 
because of the symmetry of the cluster: 9^ = 9 , for all 
( i , n) and = ±e, where the + sign applies to y, — 1, 2, 3 
and the — sign to p, = 4, 5,6. Particle-hole symmetry 
forces i) il j to vanish; however, this will no longer be the 
case in the antiferromagnetic phase. 

We probe the antiferromagnetic phase by adding to the 
cluster Hamiltonian H' the term: 


H' m = M 


EKt - n i,l) - - n i,l) 

jeA ieB 


( 20 ) 


where ni t<T = c\ a Ci t<T , A and B stand for the two sublat¬ 
tices of the honeycomb lattice, and M, the antiferromag¬ 
netic Weiss field, is an additional variational parameter. 
The values of this Weiss held on the two clusters will 
be opposite. In addition, as mentioned above, a spin- 
dependent hybridization for all (i,/i), will be 

allowed, which makes a total of 4 independent variational 
parameters in that phase. 

The matrix V(k, w) of Eq. (18) contains all the in¬ 
formation about the dispersion relation on the honey¬ 
comb lattice, including the hopping terms between the 
two clusters forming the repeated unit, as well as the hy¬ 
bridization functions associated with the baths connected 
to the two clusters. Because of the two clusters in the 
unit cell, all matrices (G', V, etc.) have a block struc¬ 
ture. The cluster Green function G' is block diagonal, 
but V isn’t. However, all the frequency dependence of 
V lies in the block-diagonal components, and all the mo¬ 
mentum dependence lies in the block off-diagonal com¬ 
ponents (that is not a general statement, but true for the 
system under study). 
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C. Limitations of CDIA and cluster methods 

The strength of CDIA, and of other cluster methods 
like VCA, CDMFT and DCA, resides in their inclusion of 
short-range spatial correlations and of dynamical correla¬ 
tions. However, they have the following limitations: (1) 
They do not take into account long range, two-particle 
fluctuations. Therefore they are insensitive to a possible 
destabilization of order by collective excitations, and in 
particular do not contain the physics behind the Mermin- 
Wagner theorem. (2) Like mean-field theory, they can¬ 
not find orderings that are not programmed into them. 
Specifically, the bath parameters or Weiss fields must al¬ 
low for a given broken symmetry to occur in order for 
the corresponding order to possibly emerge. (3) The or¬ 
der probed must be commensurate with the repeated unit 
(unit cell) of the system; incommensurate order and order 
with large periods cannot be decribed in this framework. 
(4) Anything about two-particle excitations is confined 
to the cluster itself, and suffers from strong finite-size 
effects. 

Thus CDIA will make statements about the Mott tran¬ 
sition or static (e.g. magnetic) order, but not about the 
type of spin liquid associated with the Mott phase. For 
that, other techniques must be used, as done in Ref. 17. 
Even if we were to compute the dynamical spin suscep¬ 
tibility Xij(w), it would be confined to each cluster, and 
would show a sizeable spin gap due to finite-size effects 
alone, which would lead us to the (wrong) conclusion 
that the spin liquid associated with the Mott insulator is 
short-ranged instead of algebraic. 

Both CDMFT and CDIA introduce bath orbitals to 
better capture quantum fluctuations in the time domain. 
An advantage of CDIA over CDMFT lies in the possibil¬ 
ity of adding Weiss fields to the cluster in order to probe 
broken symmetry phases more easily. Other advantages, 
described for instance in Refs 3, 23, and 26, include a 
better description of the Mott transition (with a clear 
hysteresis between the metallic and insulating solutions). 
In addition, CDMFT with a finite bath has an ambiguity 
in its self-consistency procedure, which does not exist in 
CDIA, since the latter is based on an exact variational 
principle. However, for a given cluster and bath, CDIA is 
more demanding numerically, essentially because apply¬ 
ing the variational principle requires a longer sequence of 
exact diagonalizations than the CDMFT self-consistent 
procedure, and convergence is more delicate. 


IV. RESULTS AND DISCUSSION 

We first used the CDIA to locate the Mott transition 
as U is increased from zero, for values of t' in the inter¬ 
val [ 0 , 1 ] (negative values of t' are equivalent to positive 
values, except for the chirality of Dirac points, which is 
reversed). We forbid the antiferromagnetic solution by 
setting the corresponding Weiss field to zero. The Mott 
transition is discontinuous, displays hysteresis, and oc- 
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FIG. 7. (color online) Density of states N(u) for different 
solutions on the phase diagram of Fig. 5. Top panel: t' = 0.5; 
bottom panel: t! = 0.9. Note that the two semi-metallic 
solutions, ( t ' = 0.5t, U = t) and (t' = 0.91, U = t), should 
have a vanishing density of states at the Fermi level, but this 
is hidden here by the use of a Lorenzian broadening 77 = O.Olf. 
The semi-metal with eight Dirac points (top panel) has a 
narrower gap-like feature near the Fermi level, compared to 
the semi-metal with two Dirac points only (bottom panel). 


curs at a critical value U c at which the grand potential 
B of the semi-metallic is the same as that of the insu¬ 
lating solutions. This transition line is shown as a green 
(dashed) line on Fig. 5. When comparing with previous 
results on the Mott transition in graphene , 3 one must re¬ 
call the factor of ^ appearing in front of t in the Hamil¬ 
tonian ( 1 ), which means that the scale of the U axis on 
Fig. 5 must be multiplied by 2. 

The semi-metallic side of the Mott transition is made 
of two different phases, depending on the number of dis¬ 
tinct Dirac points (2 or 8 ). This point was overlooked in 
Ref. 17. At U = 0, the transition between the two SM 
occurs at t' = l/v^, as explained above. For U > 0, 
the transition is still visible by carefully looking at the 
spectral function computed from the CDIA solution and 
is indicated by the blue squares on Fig. 5. Along the 
Mott line, this transition occurs towards t! ss 0.82. Note 
that the graphene limit (t' = 0 , indicated by a green full 
line on the figure) is singular, since the additional Dirac 
points, as well as the Berry curvature, appear as soon as 
t' 7 ^ 0 . 

We then relax the constraint on the AF Weiss field 
M to allow for long-range AF order and find the AF 
transition line shown as red squares on Fig. 5. For 
t' < 0.82f, the AF transition preempts the Mott tran¬ 
sition and the spin liquid state does not exist. But since 
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the spin-dependent hopping frustrates Neel order, the 
critical value Uaf recedes towards higher values as t' in¬ 
creases, revealing the underlying spin liquid phase when 
t' > 0.82 1. 



(c) U = 0, t' = 0.5t 




FIG. 8. (color online) Spectral function A( k, 0) at the Fermi 
level for different semi-metallic solutions in the phase diagram 
(Fig. 5). Top left: the graphene limit ( t ' = 0) at U = t, 
with two distinct Dirac points at K and K'. Top right: the 
same, in the semi-metallic solution bordering on the ASL. On 
the bottom panel, the presence of 6 additional distinct Dirac 
points is clearly visible at t' = 0.5, for both non-interacting 
(U = 0) and interacting (U = 2.5 1) solutions. The Lorenzian 
broadening is rj = 0.041 for panel (a) and r/ = O.Olt for the 
others. 


Fig. 6 shows the behavior of the Neel order parame¬ 
ter as a function of U for different values of the spin- 
dependent hopping t'. For t' < 0.4f the order pa¬ 
rameter behaves as a square root (critical exponent of 
/? = 1/2) around the critical Coulomb repulsion U c . 
This mean-field behavior occurs because cluster meth¬ 
ods do not capture long wavelength fluctuations needed 
to correctly predict critical exponents. When t' increases 
(0.4t < t' < 0.82t), the transition becomes more abrupt 
and the square root behavior disappears; we call this a 
weakly discontinuous transition. If t' > 0.82i, we observe 
clearly a jump of the order parameter which is a signa¬ 
ture of a discontinuous phase transition between the AF 
and Mott (ASL) phases. The discontinuous character 
of the transition is also seen when looking at the bath 
parameters, which show a clear jump at the transition. 

Along the boundary between the semi-metal and the 
AF phase, we observe that both the graphene Dirac 
points and the new Dirac points disappear at once. As 


soon as one enters the antiferromagnetic phase, the Weiss 
field M is nonzero. This parameter breaks parity, and it 
is precisely that symmetry that protects the graphene 
Dirac points. If the Weiss field M were part of a nonin¬ 
teracting Hamiltonian, then all Dirac point would disap¬ 
pear simultaneously, as shown around Eq. 8. Thus it is 
not unnatural for the two types of Dirac points to dis¬ 
appear together. Here the AF gap created at the Dirac 
points is a self-energy effect, but we should not be sur¬ 
prised that it affects all Dirac points simultaneously, in 
view of the U = 0 behavior when parity is broken. In 
addition, particle-hole symmetry is broken as soon as we 
enter the AFM phase. Although we are lucky enough to 
be positioned within the AFM gap by chosing fi = U/2, 
the spectral function is no longer symmetric around the 
Fermi level. 

It is remarkable that the intersection of the Neel curve 
with the Mott curve coincides with the topological tran¬ 
sition between the two types of semi-metals, even though 
the numerical procedures to determine that point are dif¬ 
ferent. Thus the four phases identified in this work meet 
at {U,t') « (4.8,0.82). 

We will not argue here why the Mott phase found at 
t' > 0.82 is an algebraic spin liquid. The argument can¬ 
not be made using CDIA results, and can be found in 
Ref. 17 and the associated supplementary material. Note 
however that this phase appears at stronger Coulomb re¬ 
pulsion ( U/t > 4.8) than in Ref. 17, because of our use of 
CDIA instead of the simpler Cluster Perturbation The¬ 
ory. 

Fig. 7 illustrates the density of states for AF insulator, 
ASL and SM regions of the phase diagram (see Fig. 5). 
The gap in the ASL and AF Mott insulator phase can 
be seen clearly at the Fermi level. For the same value 
U = t, the two semi-metallic phases have different low- 
energy structures, and each transits to a different phase 
upon increasing U. 

In Fig. 8, we represent the spectral function A(k, uj = 
0) at the Fermi level for different semi-metallic solutions 
of the phase diagram (Fig. 5). The non-interacting case 
solved analytically at the beginning is correctly repre¬ 
sented with the six additional Dirac points when t' < t' c 
as shown in Fig. 8 (c). Above this critical value and 
at t' = 0, the system is graphene-like with two Dirac 
points (Fig. 8 (a) and (b)). The features observed in (d) 
are similar to those of the noninteracting case (c), but 
broader. 

In principle, other magnetic orders could exist, and 
compete with both the simple AF order studied here and 
the spin liquid. We have probed a stripe-like collinear 
magnetic order with ordering wavevector Q = M (see 
Fig. 9), without finding any solution. Although this does 
not eliminate completely the possibility of competing or¬ 
der being present, it strengthens our point that the Mott 
phase exists when t' is close to 1 and U large enough. 







FIG. 9. (color online) Another possible magnetic order, with 
ordering wavevector Q = M, probed in this work but for 
which no solution was found in the range of U and t' covered. 


V. CONCLUSION 


We have investigated the phase diagram of the half- 
filled Kitaev-Hubbard model. The analytic solution in 
the non-interacting limit reveals a Lifshitz transition be¬ 
tween two semi-metallic states, with two and eight Dirac 
points, respectively. The Chern number is the same 


for these two phases, and the transition between the 
two occurs as the new Dirac points annihilate pairwise, 
forming a semi-Dirac point precisely at the transition . 20 
These two phases survive in the presence of interaction, 
as shown by an approximate solution of the interact¬ 
ing model using the Cluster Dynamical Impurity Ap¬ 
proximation (CDIA). In principle, the transition between 
the two semi-metals could be observed in Bloch-Zener 
tunneling 19 in a cold-atoms realization of the model. 
Overall, the phase diagram contains four phases that 
meet at a single point. On the strong coupling side, these 
are an antiferromagnetic phase at low t ', and a spin liquid 
phase (shown to be an algebraic spin liquid in Ref. 17) at 
high t'. The transition from the antiferromagnet to the 
spin liquid is discontinuous, whereas the transition from 
the semi-metal to the antiferromagnet, which pre-empts 
a Mott transition, is continuous. 

We gratefully acknowledge discussions with R. Shankar 
and M.S. Laad. Computational resources were provided 
by Compute Canada and Calcul Quebec. 
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